Test of classical nucleation theory on deeply supercooled high-pressure simulated silica 
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We test classical nucleation theory (CNT) in the case of simulations of deeply supercooled, high 
density liquid silica, as modelled by the BKS potential. We find that at density p = 4.38 g/cm 3 , 
spontaneous nucleation of crystalline stishovite occurs in conventional molecular dynamics simu- 
lations at temperature T = 3000 K, and we evaluate the nucleation rate J directly at this T via 
"brute force" sampling of nucleation events in numerous independent runs. We then use parallel, 
constrained Monte Carlo simulations to evaluate AG(n), the free energy to form a crystalline em- 
bryo containing n silicon atoms, at T = 3000, 3100, 3200 and 3300 K. By comparing the form of 
AG(n) to CNT, we test the ability of CNT to reproduce the observed behavior as we approach the 
regime where spontaneous nucleation occurs on simulation time scales. We find that the prediction 
of CNT for the n-dependence of AG(n) fits reasonably well to the data at all T studied. A/z, the 
chemical potential difference between bulk liquid and stishovite, is evaluated as a fit parameter in 
our analysis of the form of AG(n). Compared to directly determined values of A/i extracted from 
previous work, the fitted values agree only at T = 3300 K; at lower T the fitted values increasingly 
overestimate A/i as T decreases. We find that n* , the size of the critical nucleus, is approximately 
10 silicon atoms at T = 3300 K. At 3000 K, n* decreases to approximately 3, and at such small 
sizes methodological challenges arise in the evaluation of AG(n) when using standard techniques; 
indeed even the thermodynamic stability of the supercooled liquid comes into question under these 
conditions. We therefore present a modified approach that permits an estimation of AG(n) at 
3000 K. Finally, we directly evaluate at T — 3000 K the kinetic prefactors in the CNT expression 
for J, and find physically reasonable values; e.g. the diffusion length that Si atoms must travel in 
order to move from the liquid to the crystal embryo is approximately 0.2 nm. We are thereby able 
to compare the results for J at 3000 K obtained both directly and based on CNT, and find that 
they agree within an order of magnitude. In sum, our work quantifies how certain predictions of 
CNT (e.g. for Afi) break down in this deeply supercooled limit, while others [the n-dependence of 
AG(n)] are not as adversely affected. 



I. INTRODUCTION 

In recent years, computer simulations have increasingly 
been used to study the nucleation and growth of crystals 
from the supercooled liquid state. Molecular dynamics 
(MD) simulations have been particularly useful in test- 
ing, on a molecular level, the predictions of classical nu- 
cleation theory (CNT) U S S 0> la, 01 ■ In large mea- 
sure, this has been made possible by the development of 
novel computational techniques that permit the determi- 
nation, from simulations, of free energy barriers, kinetic 
prefactors, and the order parameters required to quan- 
titatively test the predictions of CNT fA H H Hoi Hl|. 
A key feature of these techniques is that they allow the 
study of nucleation under thermodynamic conditions at 
which spontaneous crystal nucleation does not occur on 
the short physical time scales accessible to conventional 
MD simulations. As a consequence, much previous work 
has focussed on testing CNT, and calculating the nul- 
cleation rate, at low to intermediate degrees of super- 
cooling, where CNT is expected to best apply. 

At the same time, spontaneous crystal nucleation is 
observed in a number of simulated liquid systems under 
very highly supercooled conditions where the nucleation 
time is comparable to or less than the simulation time 
scale (e.g. Refs. [H [H Q E EH). Current conven- 
tional MD simulations typically are able to study systems 



of a few thousand molecules over a time scale of tens of 
nanoseconds. Within these restrictions, a spontaneously 
crystallizing system will exhibit quite small crystal nuclei 
compared to those found at higher temperature T, and it 
is generally expected that CNT will not predict well the 
behavior of the system in this regime. Consequently, rel- 
atively few studies examine this deeply supercooled limit 
of nucleation behavior in the context of CNT. 

The purpose of the present work is to explore this 
deeply supercooled limit of nucleation behavior, with the 
goal of testing the limits of CNT and quantifying how 
the theory begins to fail in this regime; and also to de- 
termine the technical limits of applicability of the sim- 
ulation methods usually employed at higher T. We are 
interested in determining if it is possible to compare a nu- 
cleation rate calculated using CNT, and a rate found di- 
rectly from a spontaneously crystallizing MD simulation. 
The latter question is particularly interesting, since only 
a few simulation studies compare nucleation rates found 
from CNT to a rate calculated independently H 0, 0] , 
yet such comparisons are a key tool for developing and 
testing improved theoretical descriptions of nucleation. 

To achieve these goals, we study liquid silica as mod- 
elled by the BKS potential E3 • The thermodynamic and 
transport properties of the supercooled li quid state of this 
model have been characterized in detail |19l l2Cj . Previ- 
ous work has also evaluated the phase diagram of the 



2 



system, providing the coexistence conditions demarcat- 
ing the liquid, and several crystalline phases [2l|. Most 
significant for the current purpose, we find that the liq- 
uid spontaneously crystallizes to stishovite (2^ in our 
simulations when cooled to approximately T = 3000 K 
at density p — 4.38 g/cm 3 . The liquid at this T ex- 
hibits the two-step relaxation in its dynamical quantities 
characteristic of a deeply supercooled fluid, but it is still 
diffusive enough to reach metastable equilibrium on a 
time scale much shorter than the time scale for crystal 
nucleation. Consequently, we are able to make a direct 
calculation of the rate at 3000 K using an ensemble of 
independent MD simulations, while at the same time, we 
can determine the properties of the metastable liquid. 

We also use constrained Monte Carlo simulations of the 
liquid to calculate the free energy barrier to nucleation 
at the same density, over a range of temperatures from 
3000 K to 3300 K, to test the degree to which the predic- 
tions of CNT are satisfied on approach to T — 3000 K. 
The key predictions of CNT we wish to test relate to 
the central quantity of the theory, N(n), the equilibrium 
cluster size distribution, or the number of clusters con- 
taining n particles [g- In this work, we will track Si 
atoms only, and assume from stoichiometry that a clus- 
ter nominally of size n (n Si atoms) actually contains 3n 
atoms (n Si atoms and 2n O atoms). N(n) is inter- 
preted to yield the work AG(n) of forming a cluster of 
size n from the surrounding metastable liquid via, 



AG(n) 



In 



N(n) 



N(0) 



(1) 



where N(0) is the number of liquid-like Si atoms; so de- 
fined AG(0) = 0. Whether the distribution of cluster 
sizes is extensive or intensive (i.e., normalized or not), 
the barrier is system size independent. Within the CNT 
framework, the phenomenological model for the work is 
given by, 
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where Ap = /i s ti S h — A*Hq i s the difference in chemical 
potential between the bulk stable and metastable phases 
and a is a surface term that is proportional to the surface 
tension 7 and depends on the shape of the nuclei. At a 
critical cluster size n*, AG(n) has a maximum and clus- 
ters larger than n* will grow spontaneously, forming the 
new phase. AG(n*) then represents the free energy bar- 
rier to nucleation. In this study, we use computer simula- 
tion techniques that connect AG(n) with the probability 
of appearance of an n-sized cluster within the simula- 
tion, where the cluster is identified by a specific cluster 
criterion 0, fToL l23j | . We can then compare our barrier 
calculations with the general form suggested by Eq. [3 

According to CNT, the rate of nucleation, i.e. the rate 
at which critical nuclei go over the barrier, is 



J CNT = Kexp 



AG(r 



k B T 
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where the kinetic prefactor is given by, 
K = 24: Pn ZDn* 2/3 /\ 2 

1 crit ' 



Pn Z f crit 



(4) 
(5) 



where Z = y/\Afi\ /Qnk-QTn* is the Zeldovich factor, D 
the diffusion constant, k B is the Boltzmann constant, p n 
is the number density of particles, A is a typical distance 
particles must diffuse in order to go from the metastable 
liquid to the embryonic cluster, and /^ it is the rate at 
which particles are added to the critical nucleus. We 
note that the use of /^ it is an innovation introduced in 
Ref . 0| • In the case of diffusive barrier crossing /+ it can 
be calculated from simulation via, 



+ 

crit 



n [n*(t)-n*(0)]' 
2 



t 



(6) 



where (.} denotes an ensemble average. 

Our Monte Carlo simulations of liquid silica between 
3300 and 3000 K show that CNT describes the liquid well 
at the highest T, but that deviations in the observed and 
predicted behavior emerge at lower T. At the lowest T, 
we also identify technical difficulties associated with ob- 
taining N(n) and we describe an alternative strategy that 
at least partially addresses them. Notwithstanding these 
challenges, at the lowest T — 3000 K, we are still able to 
calculate the kinetic prefactors for the nucleation rate as 
described in CNT, so that we can compare the predicted 
rate to that calculated from direct MD simulations. De- 
spite the worsening correspondence between our results 
and the thermodynamic aspects of CNT at low T, the 
rates compare reasonably well. Whether the correspon- 
dence of the rates at this large degree of supercooling 
is peculiar to our system, or whether this is a general 
result is an open question. We also find that N(n), as 
obtained for the equilibrium system (i.e., the system that 
samples the equilibrium distribution of embryos, includ- 
ing those embryos near to, at, and beyond the nucleation 
barrier), is different from the analogous quantity for the 
metastable liquid state (i.e., the metastable equilibrium 
sampled in a conventional MD simulation prior to the 
onset of nucleation), raising the question of which distri- 
bution is more significant in determining the rate. 

In Section II, we desscribe the model system. Section 
III describes the direct MD nucleation rate calculation, 
while Section IV describes the CNT calculations and ex- 
plores the use of the metastable liquid in determining the 
free energy barrier. We subsequently present our Discus- 
sion and Conclusions. Appendix I describes our methods 
and criteria for defining a crystalline cluster. 



II. SYSTEM OF STUDY 

We study a system of 444 Si and 888 O ions governed 
by a modified BKS potential and in a cubic simulation 
cell with periodic boundaries. The modification includes 



3 



03 
Q. 

CD 



50 
45 
40 
35 
30 
25 
20 
15 
10 
5 




(a) 



1000 



5000 



4500 



4000 



' 3500 - 



3000 



2500- 



2000 



2000 



3000 4000 

T(K) 



5000 



6000 




V (cm 3 mol" ) 

FIG. 1: Location of studied liquid state points (filled circles) 
in the (a) PT and (b) VT phase diagram of BKS silica. At 
a given P and T, stishovite has a higher density than that of 
the liquid. Therefore, we also plot the stishovite state points 
corresponding to the liquid P (open squares) in (b). The 
phases shown in the diagrams are the liquid (L), /3-quartz 
(Q), coesite (C) and stishovite (S). (b) At fixed V, thermo- 
dynamic ground states are often mixtures of two coexisting 
phases. The liquid state points studied fall within the one- 
phase stability field of stishovite. Despite the proximity of our 
chosen state points to the stishovite to stishovite-plus-coesite 
boundary in the VT plane, the location in the PT plane is 
deep within the stishovite region. Dashed lines are metastable 
extensions of stishovite and /3-quartz transitions, and dotted 
lines represent the uncertainty in the location of the melting 
lines. V is given per mol ion. 



a short range additive patch to prevent "fusion" events, 
and a tapering of the real space part of the potential 
with a polynomial tail so that it smoothly reaches zero 
at 1 nm. Long range forces are handled via the Ewald 
summation. The details of the potential are given in 
Ref. IH. 




FIG. 2: Stishovite at 3000 K, as viewed down the crystallo- 
graphic c-axis. Only Si atoms are shown. 



The phase dia gram of BKS silica has been recently 
evaluated in Ref. [21j (Fig. ^) . In that work, the stabil- 
ity fields of the liquid, stishovite, /3-quartz and coesite 
have been determined. In this study we focus on the 
molar volume V = V = 4.5733 cm 3 /mol ( p = 4.3793 g 
cm~ 3 ) liquid isochore, which as shown in Fig.QIb), falls 
in the one-phase stability field of stishovite. A view of 
stishovite along the c-axis showing Si atoms only is pro- 
vided in Fig. |3 In stishovite, there are six O atoms 
surrounding each Si atom in an octahedral arrangement. 
These octahedra, connected along edges and at corners, 
arrange themselves in a compact manner. 

We perform our simulations in the NVT ensemble, 
where N is the number of molecules. Usually, nucle- 
ation experiments and simulations are done at constant 
pressure, P. However, we find that for the state points 
of interest, the critical nuclei we observe are small, and 
do not noticeably change the P of the system. Only once 
the crystallization process advances well into the growth 
stage does P or the potential energy U of the system 
change significantly. Therefore, the critical nucleus forms 
within a liquid that is characterized to a good approxima- 
tion either by the system's P or V. However, the density 
of the nucleus itself is not known. Hence, we also show 
in Fig. IHb) the V of bulk stishovite at the pressure at 
which the liquid is studied. 

P and U are needed along the Vb isochore to deter- 
mine A/x, which we obtain by extending the calculations 
described in Ref. |2l|. The diffusion coefficient is also 
needed in order to calculate the kinetic prefactor. To 
obtain these quantities, we perform MD simulations at 
constant V in both the liquid and stishovite, near and 



along Vq. First we equilibrate the system near the de- 
sired T with simple velocity scaling every 100 timesteps, 
and then we allow the system to continue in the NVE 
ensemble for about 1 ns (E is the total energy). For the 
T = 3000 K and T — 3100 K cases, numerous indepen- 
dent runs are performed, and only those that do not show 
any signs of crystallizing are used to determine desired 
quantities. 

For the liquid at T = T = 3000 K and V , the Si 
diffusion coefficient is D — 8.04 ± 0.2 x 10~ 7 cm 2 s _1 as 
determined from the slope of the mean squared displace- 
ment of Si atoms as a function of time t, and the pressure 
is found to be P = 44.0 GPa. This and higher T state 
points are shown to be in the stishovite stability field in 
the PT phase diagram [Fig.^a)]. 

At Pq, stishovite has a lower molar volume (V = 
4.364 cm 3 moD 1 ) than the melt. The stishovite state 
points corresponding to the pressure of the liquid [shown 
in Fig. ^b)] are the state points used to calculate A/z: 
A/i is calculated between the two phases at the same P, 
for which (in general) the volumes are different. 



III. NUCLEATION RATE FROM MD 
SIMULATIONS 

We use MD simulations in the NVT ensemble to cal- 
culate the nucleation rate at T and Vo in a "brute force" 
way. We first equilibrate the liquid at 5000 K, and then 
quench 198 independent configurations to 3000 K, em- 
ploying the Berendsen thermostat [24[ with a time con- 
stant of 1 ps. Other simulation details are the same as 
given in the previous section. The simulation continues 
at 3000 K until the system crystallizes. 

Fig. Ola) shows two sample time series, illustrating a 
large drop in potential energy associated with the phase 
change. Fig. E^b) shows components of the radial disti- 
bution function g(r) = gsm(r) + 2sr Si0 (r) + goo(r) for 
the metastable liquid (i.e., when the time series is stable), 
the crystallized system, and pure stishovite at T . The 
comparison of each g(r) shows that we indeed crystallize 
to stishovite. 

Fig- Et c ) shows FsiSi(q,t), the dynamic structure fac- 
tor at fixed wavenumber q obtained by considering only 
Si atoms (which diffuse more slowly than O), obtained 
from the metastable liquid portion of a simulation before 
the onset of crystallization. As a reference for the three 
wavenumbers chosen, we plot the static structure factor 
S(q) in the inset. We do not observe any time evolution 
of S(q) during the steady state liquid portions of the time 
series. With regards to the NVE simulations, we do not 
observe any significant differences in F(q,t) or S(q). 

From FsiSi(q,t), we see that the a-relaxation time for 
the system at Vq and T is approximately 100 ps. Thus, 
even though the system exhibits two-step (glassy) re- 
laxation, the relaxation time is typically much shorter 
than the nucleation times and we are able to achieve a 
metastable liquid state. 
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FIG. 3: Crystatizing liquid, (a) Determination of t x from 
the potential energy as a function of time. At t = 0, the 
thermostat is reset from 5000 K to 3000 K. When the system 
reaches a potential energy of U x = —1.81 (MJ/mol), it is 
well underway to crystallizing, (b) Structure at 3000 K as 
measured by gsisi(f) and gsio(r) (inset), of the liquid (Liq), 
the system after crystallization (X-Liq) and stishovite (Stish). 
(c) FsiSi(q,t), the dynamic structure factor for Si atoms for 
three q. Inset shows the static structure factor S(q) along with 
its components: total (circle), SiSi (square), OO (triangle) , 
and SiO (diamond). 
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If the fraction R(t) of unnucleated systems obeys a 
simple first-order rate law, then the rate of nucleation J 
can be obtained from 

In [fl(t)] = -JV(t-t Q ) , (7) 

where V is the volume of the system and to is the lag time, 
i.e. the time required to achieve a steady state of pre- 
critical nuclei [25| . However, to take advantage of Eq. [7| 
we must be able to identify the time when a particular 
system from our ensemble of runs has nucleated. Vari- 
ous methods can be used to detect crystallization such as 
examing the Voronoi volumes or counting the number of 
particle neighbors [25|. In this study, we contrast three 
approaches: first, we employ a simple energy criterion 
specific to the state point studied so that the crystalliza- 
tion time t x for a MD run is the time at which the poten- 
tial energy first reaches a value of U x = —1.81 MJ/mol. 
That very few simulation runs reached U x and then re- 
turned to a steady state liquid means that the system has 
progressed well past nucleation. In a sense, this mimics 
experimental measures which only identify a nucleation 
event by observing post critical clusters that are growing. 
Second, we allow runs to continue 500 ps after the system 
reaches a lower potential energy Ui ow = —1.82 MJ/mol. 
Using E/iowj we can also check the sensitivity of our rate 
calculation on the chosen energy threshold. 

Our third criterion is based on identifying the critical 
nucleus. In Section IV, we find the size of the critical 
nucleus, using the cluster criteria outlined in Appendix I 
to identify an n-sized cluster, to be about 3. We can then 
define a new time, t nuc as the latest time at which the 
largest cluster in the system n max < 1, i.e., the last time 
the liquid is precritical, and compare the nucleation rates 
based on our different criteria. By choosing n max < 1, 
we err on the side of making i nuc a lower bound on the 
nucleation time. 

Thus, we have three measures of the nucleation time: 
the time t x it takes to reach a potential energy U x in- 
dicative of the beginning of crystallization; the time it 
takes to reach a low energy threshold U\ ow ; and the last 
time i nuc at which the system possesses a largest cluster 
of size 1. 

Fig. 0] shows a plot of R using these three criteria R x , 
Riow(t) and i?nuc(i)j obtained using the upper and lower 
energy criteria and the critical cluster criteria respec- 
tively, as a function of time. The slope of each of these 
functions is the same within the error suggesting that our 
rate calculation is not highly sensitive to the nucleation 
criteria, while the lag time, obtained from the intercept, 
is sensitive. Using t x , we find that the 198 quenches 
from 5000 K to 3000 K yield a shortest crystallization 
time of 0.28 ns, and longest time of 10.53 ns, with an av- 
erage time of 2.12 ns. The slope of the line of best fit is 
0.60 ±0.02 ns -1 , where we have omitted data from times 
before 1 ns in the fit, and the time lag to = 0-4 ns. Given 
that our cubic box length is 2.1627 nm, this yields a rate 
of J = 0.059 ± 0.002 nm- 3 ns- 1 , or 6 x 10 34 m^s" 1 . 

We measure t nuc in each of our 198 nucleation runs 




<0 □ - 



2 4 6 8 

t (ns) 



FIG. 4: Determination of nulceation rate from R. The func- 
tion _R X (circles) is shown here to be well described (except 
for very early times) by exponential decay with rate constant 
0.60±0.02 ns -1 , determined from the line of best fit. Data be- 
tween 1 ns and 5 ns were used to obtain the fit. The functions 
fliow(i) and R n uc(t) have the same slope to within uncertainty. 
Inset: close-up at small t. 



to a resolution of 0.005 ns and find a lower estimate 
of the lag time to ~ 0.260 ns. The average time dif- 
ference is t x — t nuc — 0.212 ns, with standard deviation 

0. 117 ns. Comparing these two criteria, we find the mean 
value of the largest cluster at t x is {n max (t x )) = 39, i.e. 
about 10% of the system, with a standard deviation of 
23. Therefore, for our system, the crystallization process 
significantly lowers the energy only when about 10% of 
the system has crystallized. We also note that in about 
5% of the runs, the U x criterion is triggered prematurely, 

1. e., a low-energy fluctuation goes below U x , but then the 
system energy remains in steady state. 

We note that calculating the rate directly from the 
slope of the plots in Fig. 0] may ignore effects due to 
transient nucleation, i.e. that R decreases smoothly at 
early times instead of remaining at 1 until to |25| , so that 
our estimate of the rate and lag times are strictly lower 
bounds. Nevertheless, the independence of our slope on 
the nucleation criteria suggests that our calculation of 
the rate is robust and any corrections would be small. 

Fig. EJa) shows a critical nucleus at 3300 K (obtained 
from constrained MC simulations described later), while 
Fig. Elb) shows a postcritical crystallite of size 23 from 
an MD simulation at To. Fig. EJc) shows a representa- 
tive configuration at the end of a crystallization simula- 
tion. These pictures provide a visual confirmation that 
the liquid does indeed crystallize to stishovite; that small 
nuclei (at least for n max > 10) resemble the bulk phase; 
and that the procedure used to define clusters is able to 
track the nucleation process. 
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FIG. 5: Top: Sample critical nucleus at 3300 K containing 10 
Si atoms. Middle: A snapshot of the growing crystal embryo 
from a dynamic crystallization simulation at 3000 K when it 
contains 23 Si atoms. Bottom: Sample end configuration of 
a crystallization simulation. 



IV. CNT CALCULATIONS 

A. Free energy barrier 

The central quantity of CNT is N(n). However, it is 
not feasible to obtain N(n) through direct simulation for 
two reasons: critical clusters are typically rare, and hence 
it is difficult to gain statistics; and after the formation of 
the critical cluster, the system irreversibly evolves toward 



the crystalline state. To overcome this, we add a bias to 
the system Hamiltonian which constrains the clusters of 
interest into existence. The new constrained Hamiltonian 
is then 



He — H] 



BKS 



(8) 



where Hbks is the unbiased Hamiltonian derived from 
our BKS potential and 



(''max) — 2 (^rnax 



no) 



(9) 



is the contraint with k and uq being constants, and where 
n max , the size of the largest cluster in the system, is an 
order parameter [lj}. The N(n) measured under the con- 
straint is then related to its value in the unconstrained 
system through the relation 



(N(n)) 



{N{n) exp[0(n max )/fc B T]} 



c 



(exp [0(n max )/fc B T]) c 



(10) 



where {.) c denotes an average in the constrained en- 
semble. In the case where a cluster of size n is rare 
N(n) — P(n max ), the probability that the largest cluster 
in the system is of size n m ax, and Eg. 1101 becomes. 



(P(n max )) 



c exp [(/>(n max ) / ksT] 
(exp [</}(n max )/k B T]) c 



(11) 



It is important to note that our cluster definition ignores 
O atoms, and only uses Si atoms. Thus, a cluster of size 
n contains n Si atoms, or n S1O2 units (see Appendix). 

Since it is easier in practice to measure P(n max ) than 
N(n), we will use -P(n ma x) interchangeably with N(n) in 
the regime where the two are shown to be equal. For- 
mally, this occurs when clusters are rare, and can be 
justified by the following. Let P n be the probability that 
there is at least one cluster of size n in the system, and 
P n (i) be the probability that there are exactly i clusters 
of size n. Then, 

Pn = P„(l) + P„(2)+P„(3)... (12) 
N(n) = P n (l) + 2P n (2) + 3P n (3).... (13) 

What we mean by a rare cluster of size r is that P r (l) 
is small, and additionally that rare cluster appearance 
is independent of what other clusters are present, i.e. 
P r (2) « P r (l) x P r (l) « 0. This immediately leads 
to P n = N(n) |ToL |23| . By extension, two rare clus- 
ters of different sizes appearing at the same time also 
occurs with vanishing probability P,. +m (l) x P r (l) ~ 
[assuming P,. +m (l) < P r for m > 0, i.e. larger clusters 
are rarer], and so a rare cluster will also be the largest 
cluster in the system. From these arguments, we obtain 
P„(l) = P n = N(n) = P(n max ), for n > r (the equal- 
ity holds up to a normalization constant that is irrele- 
vant in determining the free energy). Of course, when 
N(n) 7^ P(n max ), we measure N(n) directly. 

The basic MC scheme follows that presented in 
Ref. [2^|, where short NVE MD trajectories generate 
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new configurations that are tested against the Boltzmann 
distribution. More explicitly, we begin with a configura- 
tion C\ with largest cluster n£,L. New random velocities 
drawn from the Maxwell distribution appropriate to the 
desired T are assigned to all the particles, and at this 
point the total energy is Hbks- With these new veloci- 
ties, the system evolves along a constant NVE MD tra- 
jectory for 10 timesteps, with forces derived from -ffeKSi 
to arrive at a new configuration G2 with total energy 

With a perfect inte- 



^BKS> an< ^ l ar g es t cluster n«. 



[21 

gration scheme, 
propability p given by, 



= H, 



[i] 

BKS- 



C2 is accepted with 



p = min < 1 , cxp 



fc R T 



tt[2] ttW 



(14) 



It is important to note that the acceptance criterion for 
this hybrid MD-MC method uses the total energy (kinetic 
plus potential), rather than just the potential. 

With this hybrid MC method, it is only necessary to 
evaluate the cluster size distribution of the system at 
the end of each MD mini-trajectory. The method also 
provides a way of incorporating the Ewald sums through 
multiparticle MD moves, i.e., energy changes arising from 
single particle moves are difficult to calculate efficiently 
when there are long range forces. 

In order to facilitate equilibration, we employ parallel 
tempering over a matrix of runs having different values 
of no and T. Our te mpering scheme follows the descrip- 
tions given in Refs. jlOl W\. The particulars are as fol- 
lows. Compute nodes running in parallel decide whether 
to attempt switches of configurations with neighboring 
nodes every 10 MC steps, alternating between T-switch 
and no-switch attempts. For T-switchcs, an attempt is 
made with each neighbor with probability 0.36. For no- 
switches, an attempt is made with each neighbor with 
probability 0.19. The probabilities for accepting switches 
are given in Refs. [Toll27|. In practice, it is computation- 
ally faster to switch Hamiltonians or T between proces- 
sors, rather than configurations. 

To gather data, we set up a grid of simulations with 
several values of no for each T. Simulations are seeded 
from configurations sampled from the MD crystallization 
runs at To. To initially locate n* , we set up a grid as given 
in Table HI a), with k = 8 kJ/mol. After equilibration, 
and after roughly determining the shape of the AG(n) 
curves, we set up other grids with k = 16 kJ/mol, as 
shown in Table QJb-c). 

Our results for T = 3100, 3200 and 3300 K come 
from Table IJb), and the run time is 700 000 MC steps. 
The starting configurations are those from Table [Ja). 
AG(n) is calculated over the interval from 100 000 to 
700 000 MC steps, checking that AG(n) as calculated 
separately from the intervals 100 000 - 400 000 MC steps 
and 400 000 - 700 000 MC steps do not show appreciable 
differences. For these T, N(n) ~ -P(n max ) for n > 2 and 
so we use -P(n max ) to determine AG(n). The first part 
of AG(n) is obtained by calculating N(n) directly from 
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TABLE I: Parallel tempering grids for three sets of simula- 
tions showing no and T for each node, using a parabolic con- 
straint on n max with k = 8 kJ/mol (a), and k = 16 kJ/mol 
(b-c). Each node is allowed to communicate with its near- 
est neighbor. For example, in (b), node 10 can attempt T- 
tempering switches with node 18, and no-tempering switch 
with nodes 9 and 11. 



simulations where no is small. The simulation grid in 
Table [Jc) provides a consistency check on the results for 
T = 3100 K and T = 3200 K. The run time is 350 000 
MC steps. 

Fig. El shows pieces of N(n) and -P(n max ) obtained 
from parallel simulations for T = 3200 K and for var- 
ious no- We see from the no = 1 case that for n > 2, 
-P(^max) = N(n). Therefore, only P(n max ) need be cal- 
culated for determining AG(n) beyond n — 2. Fig.Elalso 
shows the consistency of the sampling between simula- 
tions of different no near the top of the barrier. For each 
no, a portion of N(n) is recovered up to a multiplicative 
constant, or additive constant in AG(n). The pieces are 
matched using the self-consistent histogram method [23] , 
and the resulting AG(n) curve is shown in Fig. [7| The 
curves for T = 3300 and 3100 K are produced by the 
same method. 



B. Methodological challenges at T = 3000 K 

For T = 3000 K, we encounter methodological difficul- 
ties using the parabolic constraint in our MC simulations, 
apparently due to the small size of the critical nucleus. 
At this T, as wc shall see, P(n max ) ^ N(n), and so 
we must find N(n) directly, but the parabolic constraint 
together with Eq. ED does not yield adequate statistics. 
For example, in Table Cfc), node infrequently samples 
states over the barrier, and because of the large factor 
of exp [^(n max ) / ksT] in Eq. these particular states 
dominate the resulting N(n). 

To obtain N(n) at 3000 K, we replace the parabolic 
constraint with a vertical, hard wall potential by set- 
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FIG. 6: Portions of AG(n) before shifting, based on unnor- 
malized histograms for N(n) and P(n max ). These are data 
transformed via Eg. HOI and llll from the constrained into the 
BKS silica ensemble. These distributions are for T = 3200 K, 
for processors 5, 6, 7, 8 and 9 from Table db). Legend shows 
the values of no used to constrain the system. Symbols and 
bold lines indicate portions of the data used to obtain the 
complete AG(n). For the cases of no = 1 and 3, N(n) is ob- 
tained directly. The dashed line shows P(n max ) for the no = 1 
case, illustrating that already N(n) = P(n max ) for n > 2. For 
larger values of no, P(n max ) is used to obtain AG(ra). 
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TABLE II: Parallel tempering simulation grid showing T and 
limits on n max for each node, using hard wall constraints. For 
a given node, only configurations with n 1 ^^ < n max < n^ ax 
are accepted during the MC simulation. Each node is allowed 
to communicate with its nearest neighbor. For example, node 
1 can attempt T-tempering switches with node 7, and n max - 
tempering switches with nodes and 2. 

ting upper and lower bounds on n max . Any MC move 
which violates n^ ax < n max < ™max i s rejected. Eq. EH 
is modified to simply be (N(n)) = (N(n)) c for n l max < 
n max < n^ ax , i.e., N(n) is the same in the constrained 
and unconstrained ensembles within the upper and lower 
bounds on n max . The constrained Hamiltonian in this 
case becomes, 

H c = Hbks + S(n max ; n l max , n^ ax ), (15) 

where 

S(n max ) = for n l max < n max < n^ ax (16) 

oo otherwise. (17) 



15, , , . . 1 . . . . , , , , r 




n 



FIG. 7: AG(n) obtained from N(n) after piecing together re- 
sults from parallel simulations such as those shown in Fig. |S| 
Filled symbols are for data described in Table Ub), open di- 
amonds are for hard wall constraints described in Table [TTJ 
The solid curves are fits to the form given by Eq. For 
T — 3000 K only points with n < 4 are used for the fit. For 
T = 3300, a one-parameter fit is shown: A/x in Eq. [5] is ob- 
tained from independent calculations, and only oat is left to 
fit. 



We then set up the hard wall simulation as outlined in 
Table UTI taking initial configurations from TableQJc) and 
running for 1400 000 MC steps. For each node, N(n) is 
determined for n max < n < n^ ax . Small windows in n 
are use to gather good statistics as well as to prevent nu- 
cleation in the bin closest to n = 0. Error estimates are 
taken by considering different time intervals in determin- 
ing N(n). Ideally, AG(n) for T = 3100 K should be the 
same when calculated either with hard wall or parabolic 
constraints. Any discrepancy between the curves is an- 
other measure of our uncertainty in AG(n). 

As described above, we have calculated AG(n) for 
T = 3100 K with three sets of simulations, as described 
in Tables |Jb-c) and [n] In Fig. [S] we plot as crosses the 
results from TableQJc), and see that the data show good 
consistency with those obtained from Table Hfb), deviat- 
ing only beyond n = 5, the last no in Table IJc). The 
hard wall curve also shows good consistency with the 
parabolic constraint results. 

Fig.[S]also shows that while N(n) w P(n max ) holds for 
T = 3100 K, it breaks down for T = 3000 K, necessitating 
the direct calculation of N(n) which we accomplish with 
our hard wall constraints. 

However, despite the narrow binning of n max outlined 
in Table HH all the bins at T = 3000 K eventually nucle- 
ate. Fig. shows for node in Table UTI the U and rt max 
time series [panels (a) and (b)] as well as early and late 
time distributions N(n) [panel (c)] and P{n max ) [panel 
(d)]. The U time series at first glance seems stable, but 
does show larger fluctuations at later stages than at the 
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FIG. 8: AG(n) obtained from N(n) and P(n max ) for T = 
3000 K, compared to results for T — 3100 K. Open and 
filled diamonds represent the same data as in Fig. Crosses 
show results for T = 3100 K obtained from simulations in 
Table HJc), while open circles are from hard wall constraint 
simulations described in Table ITT1 All three curves agree up 
until the critical size. The line connecting plus symbols (+) 
shows for Table ITU the approximate equality of P(n max ) and 
N(n) for n > 2. The line connecting stars (*) shows that this 
equality breaks down for T = 3000 K. 



beginning. The n max time series also shows a change 
in behavior: at early times (up to 400 000 MC steps) 
n>max = or 1 is favored, while at later times n max = 1 
or 2 is favored. The early and late time N(n) profiles 
show a significant difference, as do the P(n max ) curves. 
In fact, at time beyond 500 000 MC steps, — lnP(n max ) 
monotonically decreases. 

Since n max is an order parameter, the quantity 
— fee T In P(n max ) is a free energy. The later-time curves 
shown in panel (d), therefore, would seem to indicate that 
there is no barrier to increasing n max for the system, i.e. 
that the liquid is no longer a metastable phase. Alterna- 
tively, these data might suggest that it is no longer suffi- 
cient to describe the nucleation reaction coordinate solely 
in terms of n and that an additional parameter such as 
cluster quality may be needed [28|. For example, it could 
be the case that small clusters that are well ordered may 
be over the barrier, as indicated in the second half of the 
time series, while other less ordered clusters of the same 
size are on the liquid side of the barrier. Thus, for late 
times, this bin no longer contains the metastable liquid, 
but rather a post-critical state. The late-time AG(n) 
profiles, therefore, represent a lower bound of the CNT 
barrier for this T. The early time behavior, however, still 
shows a metastable liquid according to — In P(n max ), and 
is therefore used to estimate the CNT barrier. The curves 
for T — 3100 K show no such difficulties. 

These observations highlight the difficulty in determin- 
ing AG(n) for such small values of n* , where the quality 
of the cluster may significantly affect whether a cluster of 




FIG. 9: Breakdown of n max order parameter at T = 3000 K. 
Time series are shown for node in Table [II] for U (a) and 
"■max (b). In (a), U begins in a steady state similar to the un- 
constrained MD simulations at 3000 K (0-300 000 MC steps), 
and following period of small decrease (300 000 - 500 000 MC 
steps) enters a regime with slightly larger fluctuations and 
where lower energy states are probed more often. Panel 
(b) shows a crossover near 400 000 MC steps to a regime 
where n max = 2 is favored over n max = 0. The effect on 
AG(n) = —hxN(n) + const is shown in (c), where distri- 
butions taken from different portions of the time series are 
plotted: the early time profile is distinct from the later-time 
steady state. Panel (d) plots — In P(n max ) + const, showing it 
to be monotonically decreasing at late times, apparently in- 
dicating a barrierless regime. However, it is more likely that 
the nucleation process is simply not adequately described by 
using n max alone. 



a given size is post-critical. A possible solution is to in- 
troduce another order parameter that takes into account 
the quality of the cluster when determining free energy 
profiles. 

Our main result from Fig. [H] is that we obtain 
AG(n)/k B T = 7.86 ± 0.6 for T = 3000 K. The uncer- 
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T(K) 


\Afj,\ /k B T 


|A Mflt |/fc B T 




j/ksT nm 


3000 


3.28 


5.14 


11.23 


29 


3100 (HW) 


3.12 


3.97 


10.28 


26 


3100 


3.12 


4.36 


10.84 


24 


3200 


2.96 


3.57 


10.00 


26 


3300 


2.81 


2.87 


9.08 


23 


3300 (1-par) 


2.81 


2.81 


8.96 


23 



Table CTTl Note that a = 47r(3/(47rp„)) 2 / 3 7 for spheres. 
For silica at ambient P and near T — 1500 K, experimen- 
tal values for 7 range from 0.3 to 0.7 Jm -2 , or ^y/k-^T = 
15 to 34 nm -2 30]. For comparison, the value recently 
reported for NaCl at 800 K is 7NaCi = 80 erg cm' 2 , or 
7NaCi/fcBT = 7.2 nm~ 2 . Thus, we see that at our high 
P and T, where the liquid is simpler, i.e., does not have 
a tetrahedral network, 7 is still close in value to what it 
is at ambient P, and does not have a value closer to that 
of a simple ionic liquid. Table UTTl also shows that despite 
the breakdown in the ability of CNT to predict A/j, in 
this T range, a fit of Eq. [21 to our AG(n) data still gives 
a relatively consistent estimate of 7. 

Furthermore, it is interesting to note that while 
Anat I 'ksT changes some 80% as T decreases from 3300 
to 3000 K, 7/fc B r roughly changes by only 25%. This 
perhaps indicates that the structure and/or density of 
the critical nucleus interior undergo larger changes with 
T than surface properties. 



D. Kinetic prefactor 

The crucial quantity in the kinetic prefactor is either A, 
or /+ it from Eq s.|H and|51 Following the work of Frenkel 
and co-worders [HJ, [lj , we calculate /+ it through Eq. 
Eq. follows the assumption that the addition and de- 
tachment of particles from the near-critical crystallite is 
a diffusive process. In order to measure the deviation of 
the cluster size from the critical value, i.e., the right hand 
side of the equation, we isolate 80 clusters near the crit- 
ical size from constrained MC simulations and use them 
to seed NVE simulations lasting 150 ps with randomized 
initial velocities corresponding to T — 3000 K. We then 
use multiple time origins from each time series, where at 
each time origin the configuration has n max = n* . Ad- 
ditionally, to ensure we are measuring the properties of 
clusters of critical size, each time origin is only chosen 
when the average cluster size for the preceding 1000 fs 
is between 2 and 4. Varying the averaging time or these 
upper and lower bounds does not appreciably affect the 
results. 

We plot in Fig. HOf a) the quantity 

[rt max (<) — ti*(0)] 2 ^. The plot shows a very rapid 

early time increase to a value of about 4 (inset shows 
early time behavior) or a fluctuation in size of the 
cluster of about 2 particles. Notwithstanding the early 

time change in ^[n max (£) — n*(0)] 2 ^}, we see that the 

time series enters into a diffusive regime that is linear 
in time, with |n max (t) — n*(0)| between 2 and 3. By 
fitting a line to this section, we obtain an estimate 
of the slope m = (2.0 ± 0.2) x 10 2 ns _1 which gives 
/+ it = to/2 = (1.0 ± 0.1) x 10 2 lis" 1 . This is about 3 
times larger than the value obtained for molten NaCl at 
T = 825 K and atmospheric P of 0.033 ps" 1 [HI- 

The early time behavior of ([n max (t) — n*(0)] 2 \ is 



TABLE III: Fit parameters using Eq. [5] to describe data as 
plotted in Fig. [T] The quantity |A/i| /k^T is not a fit param- 
eter, and is determined in a way described in Ref. [2l| , within 
an error of ±0.08. The label 1-par indicates a one-parameter 
fit in which only is varied. Estimates of 7 are obtained 
from aflt, assuming a spherical nucleus. 

tainty is obtained by considering different portions of the 
(early) time series when constructing AG(n). We also 
obtain n* = 3 Si atoms (SiC-2 molecules), or ~ 9 atoms 
including O. 



C. Comparison with CNT 

Fig.0shows the full AG(n) curves for the different T. 
We see from Fig.[7|that both n* and AG(n*) decrease as 
T decreases. Furthermore, we see that n* ranges from 3 
to 10. These small values make it unlikely that periodic 
boundary conditions induce catastrophic nucleation in 
our system [2^. Moreover, for T = 3000 K we find n* = 
3, and analyzing our MD simulation runs, we confirm 
that the energy (or pressure) signature of crystallization 
occurs after nucleation occurs, and that the creation of 
a critical nucleus at this T does not detectably affect the 
pressure. 

The results in Fig. [7| allow us to compare the observed 
behavior to the form of AG(n) predicted by CNT, given 
in Eq. [21 We fit Eq. |21 to the data at each T, and deter- 
mine the constants as.% and A/ig t as fitting parameters. 
These fits are shown as solid lines in Fig. [7| and show 
that the functional form of Eq. [21 satisfactorily describes 
the data at all T. For T = 3000 K, the fit is satisfactory 
only up to the top portion of the curve. 

Next, given that the presence of the critical nucleus 
does not affect the system pressure, we calculate A/x be- 
tween stishovite and the liquid at the T and P of the 
liquid VT state point under consideration (calculated as 
the difference in Gibbs free energy per mole of Si, or 
Si02 unit), based on the results presented in Ref. |21| . 
The values are presented in Table UTTl and compared to 
the corresponding values of A^ t . We find that A/Xfi t 
compares well to A/i at T — 3300 K, but that differences 
appear at lower T, getting larger as T decreases. Thus, 
though the form of Eq.|21fits the data for AG(n) at all T, 
the ability of CNT to predict A/j, is lost for T < 3300 K. 

Assuming that we have an approximately spherical nu- 
cleus, we estimate the surface tension from the fit param- 
eter oat for our range of T to be "f/k^T rs 25 nm -2 ; see 
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t (ps) 

)l 

5 10 15 

t (ps) 



Quantity 


Value 


n* 


3 


AG(n*)/k B T 


7.86 ± 0.6 


f + 

J crit 


(1.0 ±0.1) x 10 2 lis" 1 


|A M | /k B T 


3.28 ±0.08 


Pn 


43.8929 nm~ 3 


D 


(8.0 ±0.2) x 10~ 8 nn^fs" 1 


Z 


0.241 


A 


0.2 nm 


J 


6 x 10 34 m^s- 1 


jCNT 


4.1 x 10 35 m _3 s _1 


jms 


1.6 x 10 34 m" 3 s _1 



TABLE IV: Summary of calculated quantities for T = 3000 K. 




t (ps) 

FIG. 10: Calculation of f^ rit . Panel (a) shows a plot of 
([n(t) -n*(0)] 2 ) as a function of time at 3000 K. After a 
brief early time regime, the size of the cluster shows diffusive 
behavior. The slope of the line of best fit in the linear regime 
is (2.0 ± 0.2) x 10 2 ns" 1 = 2/+ it . Inset shows early time 
behavior. Panel (b) shows n max (t) for a portion of an NVE 
simulation seeded witha cluster of size n* = 3. 



plotted in the inset of Fig. HOf a). and shows a rapid 
increase corresponding to short-time fluctuations in the 
cluster size. These rapid fluctuations are seen in 
Fig. HUT b) , where we plot a representative portion of an 
NVE simulation with a critical cluster in it. Although 
short-time fluctuations can be considerable given that 
n* = 3, the general trend shown here suggests n ma x fluc- 
tuates around n*. 

All the factors required to calculate the nucleation rate 
via Eq.|3|are summarized in Table llVl The resulting rate 
is J CNT = 4.1 x 10 35 m _3 s _1 , and given the uncertainties 
in the calculated quantities, this result should be accu- 
rate within a factor of 2. Note that we have calculated 
the rate using | A/i| as obtained from independent free en- 
ergy calculations. It could be argued that |A//gt| is the 
appropriate quantity and this introduces an additional 



factor of uncertainty of y'jAjufitj / |A^| = 1.25 

The quantity A can be obtained by solving Eqs.0Jand|21 
resulting in, 

/24Dn* 2 / 3 , . 

A = \/^+— (18) 

V Jcrit 

Using our values for D, /Z it and n* , we obtain A = 
0.20 nm. To put this in perspective, the first peak of 
the Si-Si radial distribution function for the liquid at our 
state point is ~ 0.3 nm and the width of the first neighbor 
peak in the liquid gsisi(r) is about 0.1 nm. 



E. Metastable equilibrium liquid N(n) 

We now focus on the steady-state liquid that exists in 
our MD runs prior to any crystallization, which we call 
the metastable equilibrium liquid. We ask whether wc 
can extract information about the barrier to nucleation 
from this metastable equilibrium liquid, i.e., without any 
constrained MC sampling J^l). To this end, we harvest 
configurations from our direct nucleation MD runs that 
have 500 ps < t < t x — 500 ps. This reasonably ensures 
that we have configurations only from metastable liquid 
equilibrium, or steady state liquid, and have not included 
configurations that have begun to crystallize. We mea- 
sure the distributions of cluster sizes to obtain N rns (n) 
and then define the free energy 

AG ras (n) = -k B T\nN ms (n) ± const, (19) 

where the constant is chosen to ensure AG ms (0) = 
0. One might expect the two distributions N(n) and 
N ms (n) to be the same. However, N{n) is obtained by 
allowing otherwise unstable clusters to equilibrate in the 
surrounding liquid, while N ms {n) is obtained directly 
from a dynamic simulation. Furthermore, the way we 
obtain N ms {n) reduces sampling of states near the top 
of the barrier. 
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FIG. 11: Work of cluster formation at T = 3000 K 
derived from N ms (n) distributions obtained directly from 
NVT MD simulations where no crystal nucleation has oc- 
curred. Panel (a) shows AG ms (n)/fer (filled cirles) along 
with the equilibrium AG(n) /ksT (open diamonds), the curve 
(- \Afi\ n + 9.2n 2/3 )/k B T (solid curve), and a fit of the lin- 
ear part of AG rns (n)/kBT at high n. Panel (b) shows the 
quantity [AG ms (n) + \An\n\/kaT, that according to Eq. H 
should be linear function of n 2 / 3 . The line of best fit passing 
through the origin and through the data points corresponding 
to n = 1, 2, 3 and 4, has slope 9.2. 



We plot AG ms (n) in Fig. fTTTa) (filled circles). Al- 
though the first part of the data behaves like a usual bar- 
rier, the curve at larger n is straight. A linear AG ms (n) 
implies an exponential N ms (n). It also implies that the 
work required to add a particle to the cluster is indepen- 
dent of n. Clearly, this does not fit the CNT picture, 
where increasing the size of the cluster reduces the work 
required to add a particle. 

In order to extract some information from AG ms (n), 
we calculate the quantity [AG ms (n) + | A/x| n}/k-BT, which 
should be a linear function of n 2 / 3 , and plot it in 



Fig- Illf b) . The resulting curve does indeed show a lin- 
ear dependence on n 2 / 3 at small n. A fit through the 
origin and the first 4 non-zero data points yields a slope 
of 9.2. Using this slope, we construct the CNT curve 
C ms (n) = - |A^| n + 9.2n 2 / 3 and plot it in Fig. ITlTa) 
(solid line). From this curve, we obtain AG ms (n*) = 
10.74 and n* w 6. Using these parameters, while keeping 
■/crit calculated earlier, we obtain an estimate of the nu- 
cleation rate of J ms = 1.6 x 10 34 s _1 m -3 , a value closer 
to J than that obtained from N(n) (see Table EJ- 

This close agreement with J still leaves us with ques- 
tion of why AG ms (n) ^ AG(ra). In Fig.[TTl>), we plot for 
comparison AG(n), and see that is it significantly lower 
than AG ms (n). 



V. DISCUSSION 

In this study, we take advantage of a liquid state 
point that spontaneously nucleates on a time scale long 
enough to allow the determination of the properties of the 
metastable liquid. By sampling many nucleation events, 
we obtain the rate directly as shown in Fig. 0] where 
we show that the nucleation process enters a regime of 
first-order kinetics. The criterion used to determine when 
nucleation has taken place, whether an energy criteria or 
an examination of the cluster size, does not significantly 
alter our estimate of the rate. 

With the direct rate in hand, we wish to test the predic- 
tion of CNT. First, we calculate AG(n) for a series of T 
above and including T . As it is defined in Eq.QJ AG(n) 
is only formally a (relative) free energy for the case when 
n-sized clusters are rare. In such a case, AG(n) becomes 
equivalent to — In P(n max ), which is formally a free en- 
ergy. However, N(n) is the quantity of central impor- 
tance in CNT, and the interpretation of AG(n) as a free 
energy, and AG(n*) as a free energy barrier, is valid for 
more moderate supercooling or for liquid condensation 
from the supersaturated vapor. 

In this work, we probe T low enough that a consid- 
erable (small-n) portion of N(n) is not equivalent to 
P(n max ). Indeed, at T = 3000 K, where the small value 
of n* makes all cluster sizes of interest sufficiently com- 
mon, N(n) and -P(n max ) are different everywhere. This 
does not affect the formalism of CNT, merely the identi- 
fication of AG(n), as defined in Eq. as a free energy. 

This becomes important at T = 3000 K, where it is 
difficult to keep the metastable liquid from nucleating, 
even at the smallest range of n max . At this T, the free 
energy barrier according to P(n max ) is about 1.5 k^T 
[Fig. Eld)]. The possibility of a spinodal-like loss of liq- 
uid stability to the crystal becomes a prominent possi- 
bility |33. On the other hand, a set of criteria are used 
in order to distinguish liquid and crystal states. Thus, in 
a system where nucleation is taking place in a localized 
region, most of the system will still be labeled as liquid. 
Therefore, there will be more liquid-like particles (corre- 
sponding to n — 0) than single crystal-like particles (cor- 
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responding to n = 1), so that AG(n) will always show 
a barrier. Hence, if nucleation is unavoidable because 
of some spinodal process, the usefulness of interpreting 
AG(n*) as a free energy barrier is not clear. 

With these thoughts in mind, we proceed to discuss 
the progression of AG(n) with T. At T = 3300 K, 
the identification of N(n) with P(n max ) holds very well. 
The resulting AG(n) is well described by Eq. [2] Indeed, 
even the bulk value for A/z calculated independently ac- 
curately describes the data, necessitating a fit only to 
determine the surface tension term. It is important to 
note that our approximation that the appearance of the 
critical nucleus does not affect P should be weakest at 
this highest T since the critical nucleus is largest. The 
Afj, result upholds our approximation. 

The AG(n) profiles for T = 3200 and 3100 K are sim- 
ilar to the 3300 K case, except that the bulk value of 
A/n shows increasing deviation from A/ifi t . This devia- 
tion stems from either violating the assumption of cluster 
incompressibility used in the derivation of Eq. |2] 0, |33| , 
or from the (near certain) possibility that the structure 
of the small nuclei is different from bulk stishovite, and 
therefore follows a different equation of state (the bulk 
equation of state is used calculate A/i). 

Here, a note about ensembles is in order. Since we are 
at constant V, AG(n) [or — In P(n max ) for that matter] 
represent a change in Hclmholtz free energy of the sys- 
tem. The use of A/x in Eq. |2 however, is still formally 
correct [U Moreover, as we have shown, the critical 
nucleus does not significantly alter the P in our system. 
Hence, the nucleus can be regarded as being in either a 
constant V or a constant P environment. 

At T = 3000 K, clusters are insufficiently rare to iden- 
tify N(n) with P(n max ). Therefore, we calculate N(n) 
directly using hard- wall constraints, a method giving 
the same results for 3100 K as those obtained from the 
parabolic constraint. Furthermore, we see a breakdown 
in the ability of our single order parameter n max to suf- 
ficiently characterize the critical cluster. A high quality 
cluster of size 2 can be post-critical, resulting in what 
appears to be a spinodal-like profile in Fig. Eld) for late 
times. For our current purposes, we use the portion of the 
time series which explicitly shows the metastable liquid 
state in order to calculate AG(n). In order to calculate 
-P(^max) accurately at T near 3000 K, more stringent 
measures should be taken, including perhaps using an 
extra order parameter to help characterize the critical 
state better. 

Calculating the kinetic components of the CNT expres- 
sion for the rate, though more straightforward, requires 
comment. The parameter A is usually defined as the 
distance a particle must diffuse when moving from the 
surrounding fluid to the nucleating phase, which is an 
intuitive interpretation in the case of condensation in a 
dilute gas. In the case of crystal nucleation, its meaning 
is not so clear, especially when we ask how A should be 
interpreted with respect to a cluster criteria which iden- 
tifies correlations between particle environments. Never- 



theless, by calculating D and /+ it , we find A = 0.2 nm. 
This value is physically appealing, as it is less than the 
first neighbor Si-Si distance of 0.3 nm and the distance 
between the two sub peaks of the first <?siSi(?") peak for 
stishovite (the first peak is split), is about 0.1 nm. This 
reasonable value of A is evidence that CNT provides an 
adequate description of nucleation in our system, a con- 
trast to the case of molten NaCl, where A was found to 
be unphysically large [TT| . 

Our calculation of /+ it suggests there are two time 
scales associated with the dynamics of the critical cluster. 
[n max (t) — n* (0)] 2 grows rapidly over the first 40 fs before 
reaching a plateau near 4, after which, it increases at a 
much slower rate. The diffusive growth of the cluster 
occurs slowly and measurements of f^ rit on the longer 
time scale leads to reasonable values of A. The short- 
time fluctuations most likely arise arise from particles 
close to the cluster definition thresholds, for which small 
motions result in their being included or excluded in the 
crystalline cluster. 

We emphasize that despite the difficulties at 3000 K, 
the rates J and J CNT compare quite favorably, given 
similar comparisons done previously [ll|. At higher T, 
there are no difficulties with the formalism used to cal- 
culate AG(n). Furthermore, the quantitative agreement 
between AG(n) and Eq. [21 is excellent, especially given 
the fact that A/u is calculated independently for bulk 
phases. 

Another intriguing aspect of this study is the differ- 
ence between AG(n) and AG ms (n). AG(n) is obtained 
through equilibrium simulations of the constrained sys- 
tem. The constraint allows for a rigorous determination 
on N(n), allowing post-critical states to be sampled while 
determining the barrier. AG ms (n) is determined through 
MD simulations of liquid quenched from T — 5000 K to 
To, allowing the liquid to relax for several a-relaxation 
times, and collecting data only until 500 ps before crys- 
tallization is detected through the energy. Only three 
runs have t x — t nuc > 500 ps, with the largest difference 
being 781 ps. Increasing the cutoff of the time series to 
800 ps before t x does not significantly alter AG ms (n). 

Thus, AG ms (n) is constrained in a peculiar way. The 
data are pruned to include post-critical structures, so 
long as they happen to dissolve through some fluctua- 
tion. Therefore, post-critical states are sampled in a non- 
equilibrium fashion. Also, pre-critical fluctuations that 
happen to carry the system over the barrier quickly are 
not sampled well either. Therefore, near-critical states 
are sampled less often than in equilibrium [and hence 
AG ms (n) > AG(n)]. We have not determined whether 
this sampling difference is sufficient to account for the 
difference between AG" ls (n) and AG(n). 

Given these difficulties, the agreement between J and 
the rate calculated from N ms (n), following from our 
seemingly logical procedure to obtain a barrier height 
from AG ms (n) [C ms (n*)], may be fortuitous. However, 
another possibility regarding the discrepancy between 
AG ms (n) and AG(n) is that the time scale on which 
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the N(n) evolves is much longer that the a-relaxation 
time of the liquid. We note that quantities like g(r) and 
the structure factor are constant during the time 500 ps 
< t < t x — 500 ps. However, it is possible that N(n) 
evolves more slowly. In this case, the N ms (n) that we 
measure is not expected to be the same as N(n), and is 
perhaps more physically relevant in calculating the rate. 
Perhaps this is why J ms agrees better with J than does 
J CNT . Clearly, the comparison of the metastable liq- 
uid distribution and the equilibrium distribution raises 
a number of interesting questions that warrant further 
investigation. 



VI. CONCLUSIONS 

We perform NVT MD simulations of liquid silica at 
To = 3000 K and Vq — 4.5733 cm 3 /mol, corresponding to 
Po = 44.0 GPa, and calculate the rate of homogenous nu- 
cleation to stishovite to be J = (6.0 ± 0.2) x 10 34 m _3 s -1 . 
This state point is located deep in the stishovite field in 
the PT phase diagram, and within the single phase co- 
existence region of stishovite in the VT phase diagram. 
T) is at about half the melting temperature at this Po- 

We also compare this rate to that predicted by CNT. 
The work in forming a cluster of size n [AG(n)] follows 
the form predicted by CNT (Eq. EJl for T = 3000 K, 
3100 K, 3200 K and 3300 K. At 3300 K, an indepen- 
dent calculation of A/x using bulk phase values allows 
for a successful one-parameter fit of AG(n). Assuming a 
spherical nucleus, an estimate for surface tension in this 
range of T is j/k B T w 25. At T = 3000 K, the CNT form 
only fits the data well only up to n — 4, one larger than 
the critical size. Furthermore, the usual identification 
N(n) « P(n max ) breaks down strongly at this T. 

We also find that N(n) and N ms (n) differ. In fact, 
using N ms yields a CNT result that is closer to direct 
measurement. This may indicate a subtle dependence of 
the nucleation rate on the initial T from which the system 
is quenched, via slow evolution of the cluster distribution. 

Calculating the kinetic prefactor, we obtain /+ it = 
100 ± 10 ns _1 , resulting in a calculated rate of 4.1 x 
10 35 m _3 s , with an uncertainty of a factor of 2. There- 
fore CNT overestimates the rate by an order of magni- 
tude. The average distance that a Si atom must diffuse 
in attaching itself to a crystalline cluster is A = 0.2 nm. 
This length is approximately twice the width of the first 
neighbor shell of the Si-Si radial distribution function. 
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Structure 


N b 


Q 4 


Qe 


Qs 


FCC 


12 


0.19 


0.57 


0.40 


BCC 


12 


0.08 


0.54 


0.38 


HCP 


12 


0.10 


0.48 


0.32 


sc 


6 


0.76 


0.35 


0.72 


sc 


10 


0.40 


0.02 


0.60 


LIQ 


10 


0.02 


0.03 


0.02 


X-LIQ 


6 


0.21 


0.33 


0.24 


X-LIQ 


8 


0.23 


0.30 


0.27 


X-LIQ 


10 


0.23 


0.27 


0.29 


X-LIQ 


12 


0.11 


0.22 


0.33 


ST 3000K 


6 


0.39 


0.52 


0.33 


ST 3000K 


8 


0.39 


0.48 


0.35 


ST 3000K 


10 


0.40 


0.45 


0.38 


ST 3000K 


12 


0.25 


0.38 


0.42 


ST OK 


10 


0.41 


0.51 


0.42 



TABLE V: Qi for I = 4, 6 and 8, for various structures and 
choices of iV&: face-centered cubic (FCC), body-centered cubic 
(BCC), hexagonally close-packed (HCP), simple cubic (SC), 
liquid silica at T = 3000 K (LIQ), stishovite (ST) at 3000 K, 
stishovite at K, and the structure that results when the 
liquid spontaneously crystallizes to stishovite at 3000 K (X- 
LIQ). 



VIII. APPENDIX: DETERMINING CLUSTER 
SIZE 

We only consider Si atoms in our sample when deter- 
mining crystallinity. There are three reasons for this: Si 
and O atoms have different local geometry and so the 
analysis is made easier by looking only at Si-Si structure; 
it is computationally faster to do so; and our order pa- 
rameter does not track, and therefore does not influence, 
what the O atoms are doing, allowing for greater struc- 
tural freedom during the constrained MC simulations. 

In terms of determining local geometry, the difficulty 
arises in describing the different environments of Si and 
O atoms. In the case of NaCl, both species have the 
same coordination environment in the solid, and there- 
fore can be described with one scheme. In the case of 
Si02, rather than finding a way of describing Si-O, 0-0 
and Si-Si bonds separately, we choose to account for Si 
atoms only. This is also computationally faster. We as- 
sume that the strong local stoichiomctry will persist in 
the growing clusters as well. 

In order to define a crystalline cluster forming within 
the liquid, we follow the procedure laid out in |7j]. We 
need a bond order parameter that captures crystal struc- 
ture. To begin with, we use spherical harmonics Yi m (fij), 
where fy is a unit vector pointing along a bond between 
particles i and j (and thus providing elevation and az- 
imuth angles with respect to a fixed coordinate system). 
For FCC and BCC crystals, I = 6 has been used, while 
for salt (having cubic structure), I = 4 has been used [lT| . 
The first step is to define a local quantity on the particle 



15 




FIG. 12: Determination of characterization thresholds. We 
plot the probability distribution of cy values in (a) for the 
liquid, stishovite and the spontaneously crystallized liquid at 
To and Vq. The results in each case are averages from five 
configurations. Based on this plot, we choose a value of c cut = 
0.5. In (b) we plot the probability distribution of N c values 
based on c cut . At this T, stishovite appears to always have 
ten Si-Si connections per Si. Based on this plot, we choose a 
value of Nc U t = 5, at or above which a Si ion is deemed to be 
crystal-like. 



level, 

qim(i) = Y im(hj), (20) 
J'=l 

where the sum is over the Nf, bonds of particle i. A global 
measure of the overall crystallinity can be written, 

Ei=i N b(v 

and hence a quantity that does not depend on the coor- 



dinate system is, 

\ m— — l / 

Usually, Nf, is determined via a cutoff distance near the 
first minimum in the radial distribution function. Ideally, 
FCC, BCC and HCP structures have twelve neighbors, 
while simple cubic has six. However, during a simulation, 
Nf, will fluctuate. In the present work, instead of defin- 
ing a distance cut-off, we always choose the closest ten 
silicon neighbors of a given silicon atom, i.e., Nf, = 10 al- 
ways. In stishovite, the first Si-Si neighbor shell contains 
ten atoms, although the shell is split with two neighbors 
slightly closer than the other eight. 

In Table we list the Qi values for I = 4, 6 and 
8, for various crystal structures as well as for both the 
metastable liquid and the state after the liquid has spon- 
taneously crystallized (crystal with defects) . We see that 
both Q 6 and Q$ give high values for most crystals. How- 
ever, I = 8 seem to be less sensitive to the value of Nf, 
chosen. In particular, for Nf, — 10 in the case of the sim- 
ple cubic structure, Qe fares much worse that Qg. For 
our study, we do not know what the structure of pre- 
critical nuclei of stishovite is, and thus we prefer to have 
an order parameter that is more accepting of different 
structures. Therefore, we choose I = 8 

Having selected Nf, = 10, and / = 8, we now proceed to 
determine what a crystal-like atom is, and whether two 
crystal-like atoms are part of the same cluster. Having 
defined qi m in Eq.[2U| we can form a dot product c (— 1 < 
c < 1) between two neighboring Si atoms i and j, 

8 

c y = Q8m(i)qt m U), (23) 

m=-8 

where 

98m W = nyj, (24) 

(EL-8 \lSrn(i)\ 2 ) 

and q* is the complex conjugate of q. In this way, c is 
determined for every pair of neighboring atoms. For two 
atoms with very similarly oriented bonding geometry, c 
will have a value close to unity. 

The distribution of c values is plotted in Fig. ll2f a) for 
stishovite, liquid and spontaneously crystallized liquid all 
at T and Vq. We see from the plot that very few atoms 
pairs in the liquid have a value greater than about 0.75, 
while very few atom pair in stishovite have a value less 
than 0.75. However, a value of c = 0.5 provides a better 
criterion for differentiating between the liquid, and the 
spontaneously crystallized configuration. Therefore, we 
choose a cut-off value of c cut = 0.5. A pair of neighboring 
atoms i and j that have > c cu t are considered to be 
connected by a crystal-like bond. 

To define a crystal-like atom, we say that the number 
of connections iV c an atom possesses must be greater 
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than or equal to N° ut . To determine N° ut , we plot the 
distribution of N c in Fig.^Jb) for the same cases as for 
c. We see that all atoms in the stishovite crystal have 
N c = 10, while the distribution for the liquid vanishes 
near N c = 5. From the plot, any value between 5 and 



10 would serve to distinguish the liquid from the crystal. 
We choose N£ ut = 5 to be less restrictive in our choice 
of order parameter. Beyond this, clusters are defined by 
considering connections only between crystal- like atoms. 
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